Optimization of physical schemes in WRF model on cyclone simulations over Bay of Bengal using one-way ANOVA and Tukey’s test

Evaluation of appropriate physics parameterization schemes for the Weather Research and Forecasting (WRF) model is vital for accurately forecasting tropical cyclones. Three cyclones Nargis, Titli and Fani have been chosen to investigate the combination of five cloud microphysics (MP), three cumulus convection (CC), and two planetary boundary layer (PBL) schemes of the WRF model (ver. 4.0) with ARW core with respect to track and intensity to determine an optimal combination of these physical schemes. The initial and boundary conditions for sensitivity experiments are drawn from the National Centers for Environmental Prediction (NCEP) global forecasting system (GFS) data. Simulated track and intensity of three cyclonic cases are compared with the India Meteorological Department (IMD) observations. One-way analysis of variance (ANOVA) is applied to check the significance of the data obtained from the model. Further, Tukey’s test is applied for post-hoc analysis in order to identify the cluster of treatments close to IMD observations for all three cyclones. Results are obtained through the statistical analysis; average root means square error (RMSE) of intensity throughout the cyclone period and time error at landfall with the step-by-step elimination method. Through the elimination method, the optimal scheme combination is obtained. The YSU planetary boundary layer with Kain–Fritsch cumulus convection and Ferrier microphysics scheme combination is identified as an optimal combination in this study for the forecasting of tropical cyclones over the Bay of Bengal.


Scientific Reports
| (2021) 11:24412 | https://doi.org/10.1038/s41598-021-02723-z www.nature.com/scientificreports/ of combination of schemes which are significantly close to IMD observations. From this process, the treatment combinations which are significantly different from IMD observations are eliminated and thus four cluster of treatments with respect to four variables (central pressure, maximum sustained wind, longitude and latitude) are obtained for every cyclone. These clusters of scheme combinations are then tested for the error at landfall time and average Root Mean Square Error (RMSE) of central pressure (CP) and maximum sustained wind (MSW) with respect to three cyclones and through elimination process optimal scheme is obtained.

Experimental design
In order to obtain reasonable optimal combination of schemes of WRF model, three severe cyclonic cases over Bay of Bengal (Nargis-2008; Titli-2018 and Fani-2019) are considered in this study. These cyclones are severe cyclones considered as a sample for the entire population of tropical cyclones in the Bay of Bengal. The optimization method tested for the sample can then be applied over the population. The sensitivity experiments are carried out with the combination of two PBL schemes (Yonsei University scheme and Mellor-Yamada-Janjic), three CC schemes (Kain-Fritsch, Betts-Miller-Janjic and Grell-Devenyi) and five MP schemes (Ferrier, WSM6, Thompson, Lin and Kessler) on three cyclonic cases. Thus, a total of 90 sensitivity experiments (2 PBL × 3 CC × 5 MP × 3 TCs) are performed using Weather Research and Forecasting model (WRF) with ARW core, version 4.0 ( Table 1). The detailed description of WRF model physics and dynamics are given by Skamarock et al. 19 . The initial and boundary conditions are obtained from the National Centers for Environmental Prediction (NCEP) Global Analysis and Forecasting System (GFS) data available at 0.5° × 0.5° resolution, with time varying boundary conditions updated at every 6 h interval. The model is integrated over Bay of Bengal (4.37° S-27.44° N and 68.22°-107.78° E) having Mercator map projection with 9 km horizontal resolution and 35 vertical levels. The Arakawa C-grid grid staggering is used along with third-order time integration scheme and sixth-order advection scheme in both horizontal and vertical directions. Dudhia scheme 20 and the Rapid Radiative Transfer Model (RRTM) scheme 21 are used for short wave and long-wave parameterization respectively in all the sensitivity experiments. For the surface parameterization, the Noah land-surface scheme is considered. For all three cyclones, the model is initialized three days prior to the cyclone landfall and integrated upto 120 h. i.e. initial condition for cyclones Nargis, titli and Fani are considered respectively 28april2008-12UTC, 08oct2018-12UTC and 30april2019-12UTC The observed cyclone track and intensity from the India Meteorological Department (IMD) are used to estimate the track and intensity errors. In this study, central pressure (CP; hPa) and 10-m maximum sustained wind (MSW; ms −1 ) are used to estimate the intensity of tropical cyclone. In this optimality study, the best performing combination of physics schemes have been investigated individually in relation to time at landfall, intensity (CP, MSW) and track (latitude, longitude) along the trajectory of cyclone. The physics schemes as well as their codes are presented in Table 1. The treatments are coded from 0 to 30, '0' being the IMD observation. Treatment 1-30 are combinations ordered by considering the PBL scheme at first followed by CC scheme and MP scheme viz. treatment 1 will be a combination of first PBL-YSU with first CC-KF and MP Ferrier. The treatments 1 to 5 are formed by keeping the first PBL and first CC scheme constant and orderly varying the five MP schemes listed in Table 1. The treatments 6-10 are created with YSU PBL, second CC scheme viz. BMJ and varying the MP schemes. Similarly, all the other treatment combinations are created. The means of 31 treatments are compared in order to check the variability in the average effect of treatments using one-way ANOVA. Further, groups of treatments which provide values significantly close to treatment '0' are determined using Tukey's post-hoc test, while the other treatments are eliminated. The RMSE of CP and MSW are computed every 6 h and averaged over the entire duration of cyclone. The clusters of treatments obtained from the Tukey's test are checked for RMSE and time error at landfall. The treatment which is consistently close to treatment 0 (IMD observations) for all the variables (CP, MSW, longitude and latitude) with respect to three cyclones with the elimination procedure, is termed as the optimal combination of scheme.

Results
Titli. One-way ANOVA provides results from which we can identify whether the average effect of treatments is similar to one another. For cyclone Titli, this analysis (Tables 2 and 3 Tables 4 and 5). For Titli, the group of treatments significantly close to treatment 0 for CP are 16, 9, 6, 12, 14, 1, 4, 5, 20 and 13; and MSW are treatments 10, 29, 28, 24, 22, 27, 1, 19, 2, 18, 4 and 3 respectively (Tables 14b and 15b). The treatments that are much the same as IMD observations       (Fig. 1). The time error at landfall is accurately forecasted by treatments 1 and 2 for cyclone Nargis (Fig. 2). Overall, the treatment 1 is consistent in the case of Titli, Fani and Nargis.

Discussion and conclusion
One way ANOVA has been conducted on the values obtained after the model simulation with thirty combinations of physics schemes (treatments 1-30) along with the IMD observations (treatment 0) for intensity and track variables. The results obtained from the one-way ANOVA shows that the average effect of treatments for all the variables (intensity and track) are significantly different between themselves. The observed significance level is less than α = 0.10 in the case of intensity and track variables for cyclone Titli, which suggests that average effect of all treatments can be differentiated within the considered CI. On the contrary, it is noticeable for cyclone Fani, average effect of all the treatments differs significantly for longitude, CP and MSW. The treatments for latitude variable are similar to one another in the 90% CI. The p values for CP and MSW of cyclone Nargis are less than α = 0.10; which implies that all the treatments differ significantly within themselves as well as with the observations. One way analysis of variance (ANOVA) rejects the null hypothesis for latitude whereas it accepts the same for longitude. This signifies that the treatment means are significantly different among themselves for latitude and do not vary in case of longitude. The basic purpose of the study is to identify which treatment combination is close to IMD observations. Therefore, Tukey's test is used for post-hoc analysis. Significant difference in mean effect of the treatments concludes that their simulated population means were equal at 90% confidence interval, and it was difficult to distinguish them from one another. Tukey's post-hoc test is applied individually on the same variables for average effect of all treatments. This post-hoc test uses the "Honest Significant Difference", a number that represents the distance between groups, to compare every mean with every other treatment mean. Less the difference, the more similar the treatments are to one another. In case of cyclone Titli, Tukey's test analysis for CP exhibits least differences for treatment 9 and 6 with treatment 0, followed by treatment 12, 16, 13, 14, 20, 5, 4 and 1. Whereas, the treatments 10, 29, 28, 24, 22, 27, 1, 19, 2, 18, 4 and 3 are found to be substantially close to IMD observations for MSW (Table 14b). In case of latitude, the posthoc analysis reveals that treatment 1, 4, 2, 3, 16 and 23; while for longitude the treatments 1, 12, 28, 14 and 3 are significantly close to treatment 0 respectively (Table 15b)  www.nature.com/scientificreports/ The analysis of CP for cyclone Nargis displays treatment 18 to be closest to treatment 0, followed by treatment 2, 3, 19, 4, 5 17, and 1 (Table 14c). But it does not imply that 5, 17 and 1 are less significant. For wind intensity, it may be observed that treatment 2 and 18 are closest to observations. Yet, treatments 3, 4, 9, 1 and 17 are considerably comparable to treatment 0. The treatments 20 and 5 for longitude and, treatments 9 and 8 for latitude have minimal difference with treatment 0. The group of treatments 18, 2, 3, 19, 4, 16, 25, 6, 17, 7, 10 and 1 are also in close proximity to IMD observations for longitude. Whereas, for latitude, the treatments 13, 2, 1, 23, 5 and 20 are similar to treatment 0. At this stage, the treatments which are inconsistent in all the variables are eliminated. Out of the remaining treatments, the treatments 1 and 5 for CP and treatments 1 and 4 for MSW are having average RMSE significantly close to IMD observations respectively. Accuracy in time of landfall is another vital factor to be considered while optimization. The treatments 1 and 2 have provided accuracy in the time error at landfall for cyclone Nargis. Yet, treatment 2 is eliminated as it does not comply with the RMSE in comparison to treatment 1. This makes treatment 1 the optimal scheme combination for cyclone Nargis. The present results strongly corroborate with the earlier studies on tropical cyclone simulations using the combination YSU PBL, KF CC and Ferrier MP 9,10,22 . Overall, based on the elimination process treatment 1 (YSU-KF-Ferrier) exhibits significantly close and persistence with the IMD observations.

Methods
The process of scrutinization of different combination of three physics schemes of WRF model is discussed in this section. The scheme combinations are accessed with respect to the intensity variables (CP; MSW) and track (longitude, latitude) individually on three cyclonic cases over the Bay of Bengal. Here, the tropical cyclone parameters (CP, MSW, longitude, latitude) represent the factors, and the scheme combination represent the treatments. The intensity variables are considered as they play an important role in determining the strength and stage of tropical cyclone whereas the position (longitude, latitude) is required to accurately depict the movement of cyclone. These four factors help the policy makers to predict the area of landfall, and the damage as a result of this catastrophe.
The 0th treatment is the IMD observation which is followed by the model results obtained from 30 combinations (treatments 1-30), which accounts to 31 treatments. The optimization algorithm to obtain the optimal physics scheme combination employed in the study are described as follows.
(1) Firstly, one-way analysis of variance (ANOVA) is applied on 31 treatment combinations. Then significant difference has been obtained among these treatments for each variable (CP, MSW, longitude and latitude) corresponding to each cyclone; and it is found that there is a significant difference among the 31 treatment combinations. (2) Since there was a significant difference among the 31 different treatments so, pair of treatments which are significantly different are found out through Tukey's test. This test has been applied individually on four variables of three cyclones respectively. (3) The treatments (treatments 1-30) which were consistently having no significant difference with IMD observations (treatment 0) were then tested for RMSE (CP and MSW). Average RMSE was computed over the whole TC period and the group of treatments which have significantly minimum values are found. Further, landfall time was also tested, and the treatments have been selected which were nearer to IMD observations. (4) Consistence occurrence of the treatment in all the above steps is treated as the optimal treatment. To test the null hypothesis, One-way Analysis of Variance (ANOVA) is applied for all the factors and on all three tropical cyclone cases (Nargis, Titli and Fani). ANOVA uses F-test for statistical significance, which allows for comparison of multiple means at once. Thus, the error is computed for the whole set of comparisons rather than individual. If any of the treatment means is significantly different from the overall mean, then the null hypothesis is rejected. The hypothesis is checked at 90% confidence interval (CI). The p value is then compared to the critical value p α at degrees of freedom (k, N − k), where α is the significance level i.e., 0.10; k is the number of treatments, and N is the total sample size. If the p value is less than p α (k, N − k), then the null hypothesis is rejected.
If the null hypothesis is rejected then ANOVA identifies whether there are differences among the levels of individual variables viz. differences among the treatments for individual variables, but it will not specify which differences are significant. With the purpose to identify the difference between the treatments, Tukey's Post-Hoc test is applied to test the following hypothesis.
H 0 There is no significant difference between the average effect of any two treatments. H 1 There is a significant difference between the average effect of any two treatments. Tukey's test runs pairwise comparisons among each of the groups, and uses a conservative error estimate to find the groups which are statistically different from one another. The output of the Tukey's test provides a table where treatments are arranged in order of error, from where one can identify the treatment closest to 0th viz. IMD observations. This post-hoc analysis is conducted respectively for all the four variables corresponding to the three cyclone cases.
From the results obtained from the Tukey's test, a cluster of different treatments close to observations are obtained for individual cyclonic cases with their respective variables. Treatments that are identified close to observed values need to be checked for optimality. Since, inaccuracy in forecasting of landfall time error and its intensity could result in devastation of life and infrastructure over a region thus, the RMSE computed for the CP, MSW, and the time of landfall are looked into for respective treatments. The treatment which is consistently close to IMD observations in the cyclones is suggested as the optimum treatment through elimination process for the simulation of tropical cyclones over the Bay of Bengal ("Supplementary information").
Statistical analysis. Statistical analyses are conducted using SPSS, version 22.0 (SPSS, Inc., Chicago, IL, United States). One way ANOVA is used to check the significance of observations obtained from model simulations with various combination of physics schemes hereafter termed as treatments. Results from ANOVA are reported with 90% confidence intervals and are deemed significant at p < 0.10. Tukey's post-hoc test is further used to test the significant difference between any two treatment means.